function [u_GPSR_BB_Cont,mse_BB_mono_cont]=IST_CS(Fp,Fmask,opts)

[m n] = size(Fp);
scale=sqrt(m*n);

% definde the function handles that compute 
% the compressive sensing matrix
R = @(x) fft2(x).*Fmask/scale;
RT = @(x) real(ifft2(x.*Fmask*scale));

% define the function handles that compute 
% the products by W (inverse DWT) and W' (DWT)
% other tunable parameters 
if ~isfield(opts,'Jwin') opts.Jwin=floor(log2(m))-1;  end   
if ~isfield(opts,'wav') opts.wav= daubcqf(8);  end  

wav = opts.wav;
W = @(x) midwt(x,wav,opts.Jwin);
WT = @(x) mdwt(x,wav,opts.Jwin);

%Finally define the function handles that compute 
% the products by A = RW  and A' =W'*R' 
A = @(x) R(W(x));
AT = @(x) WT(RT(x));
I0=AT(Fp);
% regularization parameter
if ~isfield(opts,'tau')  opts.tau=0.0005*max(abs(I0(:))); end %% regularization term scale
if ~isfield(opts,'first_tau_factor')  opts.tau= 0.8*(max(abs(I0(:)))/opts.tau); end %% regularization term scale

tau=opts.tau;
%first_tau_factor = 0.008*(max(abs(I0(:)))/tau);
first_tau_factor = opts.first_tau_factor;

% regularization parameter


% Run IST until the relative change in objective function is no
% larger than tolA
% [theta_ist,theta_debias,obj_IST,times_IST,debias_s,mses_IST]= ...
% 	 IST(Fp,A,opts.tau,...
% 	'Debias',1,...
% 	'AT',AT,... 
%     'True_x',WT(opts.I),...
% 	'Initialization',AT(Fp),...
% 	'StopCriterion',1,...
% 	'ToleranceA',tolA);
% 
% 
% % Now, run the GPSR functions, until they reach the same value
% % of objective function reached by IST.
% [theta1,theta_debias1,obj_GPSR_Basic,times_GPSR_Basic,debias_s,mses_GPSR_Basic]= ...
% 	GPSR_Basic(Fp,A,opts.tau,...
% 	'Debias',1,...
% 	'AT',AT,... 
%     'True_x',WT(opts.I),...
% 	'Initialization',AT(Fp),...
% 	'StopCriterion',4,...
% 	'ToleranceA',obj_IST(end));
% 
% [theta2,theta_debias2,obj_QP_BB_notmono,times_QP_BB_notmono,debias_s,mses_QP_BB_notmono]= ...
% 	GPSR_BB(Fp,A,opts.tau,...
% 	'AT', AT,...
% 	'Debias',1,...
% 	'Initialization',AT(Fp),...
%     'True_x',WT(opts.I),...
% 	'Monotone',0,...
% 	'StopCriterion',4,...
% 	'ToleranceA',obj_IST(end));
% regularization parameter
steps = 200;

debias = 1;

stopCri=3;
tolA=1.e-5;

disp('Starting GPSR BB monotonic with continuation')
[x_BB_mono_cont,x_debias_BB_mono_cont,obj_BB_mono_cont,...
    times_BB_mono_cont,debias_start_BB_mono,mse_BB_mono_cont]= ...
         GPSR_BB(Fp,A,tau,...
         'Debias',debias,...
         'AT', AT,...
         'Continuation',1,...
         'ContinuationSteps',steps,...
         'FirstTauFactor',opts.first_tau_factor,...
         'Monotone',1,...
         'Initialization',I0,...
         'True_x',WT(opts.I),...
         'StopCriterion',stopCri,...
       	 'ToleranceA',tolA,...
         'Verbose',0);
t_BB_mono_cont = times_BB_mono_cont(end);

if prod(size(x_debias_BB_mono_cont))~=0
  u_GPSR_BB_Cont=W(x_debias_BB_mono_cont);
else
 u_GPSR_BB_Cont=W(x_BB_mono_cont);
end
% % Now, run the GPSR functions, until they reach the same value
% % of objective function reached by IST.
% [theta,theta_debias,obj_GPSR_Basic,times_GPSR_Basic,debias_s,mses_GPSR_Basic]= ...
% 	GPSR_Basic(Fp,A,tau,...
% 	'Debias',1,...
% 	'AT',AT,... 
%     'True_x',WT(f),...
% 	'Initialization',AT(y),...
% 	'StopCriterion',4,...
% 	'ToleranceA',obj_IST(end));
% 
% [theta,theta_debias,obj_QP_BB_notmono,times_QP_BB_notmono,debias_s,mses_QP_BB_notmono]= ...
% 	GPSR_BB(y,A,tau,...
% 	'AT', AT,...
% 	'Debias',1,...
% 	'Initialization',AT(y),...
%     'True_x',WT(f),...
% 	'Monotone',0,...
% 	'StopCriterion',4,...
% 	'ToleranceA',obj_IST(end));
% 
% [theta,theta_debias,obj_QP_BB_mono,times_QP_BB_mono,debias_start,mses_QP_BB_mono]= ...
% 	GPSR_BB(y,A,tau,...
% 	'AT', AT,...
% 	'Debias',0,...
% 	'Initialization', AT(y),...
%     'True_x',WT(f),...
% 	'Monotone',1,...
% 	'StopCriterion',4,...
% 	'ToleranceA',obj_IST(end));
% 
% 



